\title{Matemati\v{c}no modeliranje -- 3. Doma\v{c}a naloga}
\author{Ga\v{s}per Urh\\Jon Natanael Muhovi\v{c}\\Matev\v{z} Ropret\\Miran Lipova\v{c}a}
\date{\today}
\documentclass[12pt]{article}
\usepackage{thumbpdf}
\usepackage[slovene]{babel}
\usepackage[utf8]{inputenc}
\usepackage{amsmath}
\usepackage{listings}
\usepackage[pdftex]{graphicx}
\usepackage{float}
\begin{document}
\maketitle
\thispagestyle{empty}
\newpage
\thispagestyle{empty}
\tableofcontents

\newpage

\section{Opis problema}
Na"sa naloga je bila, da za dane vrednosti $u_{i,j} := u(x_i,y_j)$ implementiramo $C^2$ interpolacijsko funkcijo $u(x,y)$ na pravokotniku $[x_1,x_m] \times [y_1,y_n]$.

\section{Na"cini interpolacije}
Obstaja mnogo na"cinov interpolacije ploskve, od teh je najenostavnejši kar linearna interpolacija
med točkami. Cilj tega poročila je predstavitev interpolacije, s katerim dobimo
gladko ploskev po interpolaciji.

\begin{figure}[H]
\centering
\includegraphics*[scale=0.65]{sincos_linearna.png}
\caption{Linearno interpolirana ploskev za enačbo $z=sin(x)+cos(y)$}
\label{fig:sincos_lineinterp}
\end{figure}


Za la"zje razumevanje bo teorija predstavljena v 2D koordinatnem sistemu.

Vzemimo tocki $A$ in $B$, ki ležita na osi $X$ (os $X$ in $Y$ tečeta od 0 do
1), ter točki $Y_A$ in $Y_B$, ki le"zita na $Y$ osi.
Točka $Y_A$ naj bo višina grafa v točki $A$, točka $Y_B$ pa višina v točki $B$. Poiskati moramo funkcijo, ki
nam višino $Y_A$ preslika v višino $Y_B$, ko se premikamo po osi $X$ iz točke $A$ proti $B$.
Najenostavnejša funkcija je, kot je omenjeno zgoraj, kar linearna interpolacija. 

Torej moramo za linearno interpolacijo poiskati tako utežno funkcijo, ki ima v $A$ višino $Y_A$, v $B$ višino $Y_B$, vmes pa se spreminja linearno v
odvisnosti od nekega parametra, ki ga poimenujmo $t$.
Vzemimo parameter $t$, ki teče od 1 do 0 (v točki $A$ ima vrednost 1, v točki $B$ pa 0). Hitro vidimo,
da utežna funkcija $u(t)=Y_A t + Y_B (1-t)$ zadosti pogoju za linerano interpolacijo.

Sedaj zgolj združimo interpolirane vrednosti v robnih točkah z njenimi sosedi in dobimo
linearni zlepek (oz linearno interpolirano površino v 3D).

Seveda če hočemo imeti gladko ploskev ta način interpolacije očitno ne zadošča.



\section{Gladka ute"zna funkcija}
\label{sec:utezna}


Če hočemo imeti res gladko ploskev, morajo imeti sosednje točke enake višine ter prve in druge odvode.
Torej moramo ugotoviti neko funkcijo, ki deluje na podobnem principu kot funkcija za linearno interpolacijo,
ampak mora zadoščati pogojem, da ima v robnih točkah enake prve in druge odvode (ti odvodi naj imajo vrednost 0, da se simetrično ujema z obeh strani).

Izkaže se, da lahko s polinomom 5. stopnje simuliramo ravno tako utež.

\begin{figure}[H]
\centering
\includegraphics*[scale=0.45]{q.png}
\caption{Modra krivulja je graf utežne funkcije $q$. $x$ teče od 0 do 1, rdeča premica pa seka ta graf natanko na polovici ($x=0.5$)}
\label{fig:qPolinonUtez}
\end{figure}


S pomočjo rdeče premice vidimo, da je graf obratno-simetričen. 
Če graf prerežemo v točki, kjer ga seka rdeča premica in ga zavrtimo za 180 stopinj vidimo,
da popolnoma prekrije drugi del grafa, kar nam vizualno pove, da je res simetričen. Ta simetričnost
na obeh straneh zagotavlja, da se bo vrednost, interpolirana s to funkcijo, lepo počasi spreminjala ob 
robnih točkah ter da ne bo ostrih robov kot pri linearni interpolaciji, saj enakost prvih in drugih odvodov
zagotavlja, da se bo višina končne točke trenutnega zlepka počasi izravnala, kar nam da vtis, da je zlepek povsem gladek ob prehodu.


Funkcijo za tako utež lahko predstavimo s polinomom 5. stopnje, ki izpolnjuje pogojem, podanim v nalogi. Torej imamo polinom naslednje oblike:

    $$ q(x) = ax^5  + bx^4 + ex^3 + cx^2 + dx + f $$

Vidimo, da imamo 6 neznank v ena"cbi, zato upo"stevamo "se vrednosti funkcije na robovih.

\begin{eqnarray*} 
   q(1) &=& 0 \\
   q(0)&=& 1\\
   q'(1) &=& 0 \\
   q'(0)&=& 0\\
   q''(1) &=& 0 \\
   q''(0)&=& 0\\
\end{eqnarray*}

 Kar pomeni, da re"sujemo sistem ena"cb s 6 neznankami in 6 ena"cbami. Za re"sitev dobimo vrednosti:

\begin{eqnarray*} 
   a &=& -6 \\
   b&=&15\\
   c&=&-1\\
   d&=&0\\
   e&=&0\\
   f&=&1
\end{eqnarray*}

Na"s kon"cni polinom je:

$$ q(x) = -6x^5 + 15x^4 - 10x^3 + 1 $$


Ko imamo enkrat koeficiente je implementacija funkcije \texttt{utez} enostavna:
\lstset{language=Matlab}

\vspace{10pt}
\scriptsize
\begin{lstlisting}
function [w]=utez(x)
        w=-6*x^5 + 15*x^4 - 10*x^3 + 0*x^2 + 0*x^1 + 1;
endfunction
\end{lstlisting}
\normalsize
\vspace{10pt}

Na sliki \ref{fig:sincos_qSamoZ} lahko vidimo, da se pri interpolaciji ploskve opazijo valovčki, ki so podobni
utežni funkciji $q$. Ta rezultat sicer ni slab, a je v mnogih primerih zaželjena še bolj gladka ploskev.

\begin{figure}[H]
\centering
\includegraphics*[scale=0.75]{sincos_qSamoZ.png}
\caption{Ploskev z enačbo $z=sin(x)+cos(y)$, katere višine so interpolirane zgolj z utežno funkcijo $q$}
\label{fig:sincos_qSamoZ}
\end{figure}



\section{Ra"cunanje odvodov za poljubno površino}
V primeru, da imamo podano naključno površino, ne moremo izračunati odvodov v znanih točkah $(x_i,y_i)$, zato jih moramo v teh to"ckah dolo"citi s približki. Za to uporabimo formule, ki so podane v navodilih. Recimo, da imamo podano površino z neznano funkcijo $u(x,y)$. Približka za prva odvoda $u_x$ in $u_y$ se izračuna na sledeči način:

\begin{eqnarray*}
u_x(x_i,y_j) &=& \frac{u_{i+1,j} - u_{i-1,j}}{2h} \\
u_y(x_i,y_j) &=& \frac{u_{i,j+1} - u_{i,j-1}}{2h}
\end{eqnarray*}

kjer je $h$ razdalja med sosednjima to"ckama v mre"zi. Drugi odvodi:

\begin{eqnarray*}
u_{xx}(x_i,y_j) &=& \frac{u_{i+1,j} - 2u_{i,j} + u_{i-1,j}}{h^2} \\
u_{yy}(x_i,y_j) &=& \frac{u_{i,j+1} - 2u_{i,j} + u_{i,j-1}}{h^2} \\
u_{xy}(x_i,y_j) &=& \frac{u_{i+1,j+1} - u_{i-1,j+1} - u_{i+1,j-1} + u_{i-1,j-1}}{h^2}
\end{eqnarray*}

Koda za računanje odvoda po $x$, dvakratnega odvoda po $x$ ter mešanega odvoda po $x$ in $y$ 
(\texttt{indexX} in \texttt{indexY} sta indeksa za v matriko višin).
\vspace{10pt}
\scriptsize
\begin{lstlisting}

  u_x=m(indexY,indexX+1)-m(indexY,indexX-1);
  u_x=u_x/(2*h);
  
  u_xx=m(indexY,indexX+1)+m(indexY,indexX-1)-2*m(indexY,indexX);
  u_xx=u_xx/(h*h);
  
  u_xy=m(indexY+1*yIndexFac,indexX+1)-m(indexY+1*yIndexFac,indexX-1) ...
       -m(indexY-1*yIndexFac,indexX+1)+m(indexY-1*yIndexFac,indexX-1);
  u_xy=u_xy/(h*h);

\end{lstlisting}
\normalsize
\vspace{10pt}

\section{Naša metoda za gladko interpolacijo površine}
Zdaj, ko lahko dobimo vrednosti funkcije in njenih odvodov na glavnih to"ckah mre"ze
se osredoto"cimo na dolo"canje interpolacijske funkcije na posameznem kvadratu.
Definirajmo funkcije $f_{00}(x,y)$, $f_{10}(x,y)$, $f_{01}(x,y)$ in $f_{11}(x,y)$, kjer $x$ in $y$ tečeta od 0 do 1.
Vsaka od teh točk je na enem od ogljišču kvadrata, katerega vmesne točke moramo gladko interpolirati.


\begin{figure}[H]
\centering
\includegraphics*[scale=0.45]{kvadrat.png}
\caption{Razpored funkcij $f_{ij}$ po kvadratu}
\label{fig:kvadrat}
\end{figure}


Funkcije $f_{ij}$ imajo obliko:
$$ f(x,y) = ax^2 + by^2 + cxy + dx + ey + f $$
Funkcija je take oblike zato, da je gladka in da jo lahko rešimo s poznanimi količinami, kot bomo videli. Druge odvode potrebujemo zato, da dosežemo popolno gladkost ploskve (če bi upoštevali samo prvi odvod bi pri senčenju ploskve opazili, da ploskev v stičiščih zlepkov ni povsem gladka).
Če to funkcijo odvajamo po $x$ ($dx$), dvakrat po $x$ ($ddx$), po $y$ ($dy$) in dvakrat po $y$ ($dyy$) ter še mešano po $x$ in $y$ ($dxy$) dobimo:
\begin{eqnarray*}
f_{x}(x,y) &=& 2xa+cy+d \\
f_{y}(x,y) &=& 2by+cx+e \\
f_{xx}(x,y) &=& 2a\\
f_{yy}(x,y) &=& 2b \\
f_{xy}(x,y) &=& c
\end{eqnarray*}
Primer: v točki (0,0) dobimo, naslednje parametre:
\begin{eqnarray*}
a &=& \frac{f_{xx}(0,0)}{2}\\
b &=& \frac{f_{yy}(0,0)}{2} \\
c &=& f_{xy}(0,0)\\
d &=& f_{x}(0,0)\\
e &=& f_{y}(0,0) \\
f &=& f(0,0)
\end{eqnarray*}


Vsaka funkcija $f_{ij}$ moramo biti definirana tako, da ima v svoji točki isto vrednost kot
višina podane točke v tistem ogljišču kvadrata. Iz tega sledi, da morajo biti vsi členi v 
tej točki 0, razen $f$, ki mora imeti isto vrednost kot je dejanska višina v tej točki.
Iz tega lahko tudi ugotovimo, kako mora izgledati funckija $f_{ij}$ za vsako posamično točko:


\begin{align*}
 f_{00}(x,y) &= ax^2 + by^2 + cxy + dx + ey + f \\
 f_{10}(x,y) &= a(1-x)^2 + by^2 + c(1-x)y + d(1-x) + ey + f \\
 f_{01}(x,y) &= ax^2 + b(1-y)^2 + cx(1-y) + dx + e(1-y) + f \\
 f_{11}(x,y) &= a(1-x)^2 + b(1-y)^2 + c(1-x)(1-y) + d(1-x) + e(1-y) + f \\
\end{align*}

Ker poznamo vrednosti $x$ in $y$ v vseh robnih točkah, lahko te funckije brez težav definiramo tako,
da jim bodo v svojem ogljišču odpadli vsi členi razen $f$.

Kot primer vzemimo funkcijo

  $$ f_{10}(x,y) = a(1-x)^2 + by^2 + c(1-x)y + d(1-x) + ey + f $$

ki leži 
v ogljišču s koordinatama $(x,y)=(1,0)$. Če v funkcijo  $f_{10}$ vsavimo $x=1$ in $y=0$ se uničijo
vsi členi razen $f$, tako da je v tej točki funkcija enaka  $ f_{10}(x,y) = f $.

Tako obnašanje funkcij $f_{ij}$ želimo, saj zaradi tega ostane višina v robnih točkah enaka.

Sedaj poznamo vrednost parametra $f$, a moramo še poiskati vrednosti parametrov $a$,$b$,$c$,$d$ in $e$.
Potrebujemo torej še 5 enačb, saj nam je ostalo še 5 neznanih parametrov.
Zgoraj smo omenili, da poznamo odvode $dx$,$ddx$,$dy$,$ddy$ in $dxy$ v vsaki točki $f_{ij}$. Torej če vsako od
funkcij $f_{ij}$ odvajamo po $x$, dvakrat po $x$, $y$, dvakrat po $y$ in $xy$ dobimo ravno še 5 enačb, katerih približne vrednosti 
pa smo že izračunali s pomočjo sosednjih točk. Sedaj lahko rešimo sistem enačb za vsako od $f_{ij}$.

Primer kode z  rešitvijo sistema za funkcijo $f_{10} $ (parameter $f$ je podan kot argument $z$):

\vspace{10pt}
\scriptsize
\begin{lstlisting}

function [v]=f10(x,y,z,xv,yv,matrikaVisin,xi,yi,velKv)
  x=1-x;
  [dx,dy,ddx,ddy,dxy]=izracunajOdvodeZaDanoTocko(matrikaVisin,xi,yi,velKv);  
  a=ddx/2;
  b=ddy/2;
  c=-dxy;
  d=-dx;
  e=dy;
  v=a*x^2+b*y^2+c*x*y+d*x+e*y+z;
endfunction


\end{lstlisting}
\normalsize
\vspace{10pt}

\section{Interpolacija}
Ko imamo enkrat ute"zi in funkcije $f_{00}$, $f_{10}$, $f_{01}$, $f_{11}$, nam "se preostane,
da implementiramo interpolacijsko funkcijo kot ute"zeno vsoto teh osnovnih funkcij oblike.

\begin{eqnarray*}
    u(x,y) = f_{00}(x,y)q(x)q(y) + f_{10}(x,y)q(1-x)q(y) \\+ f_{01}(x,y)q(x)q(1-y) + f_{11}(x,y)q(1-x)q(1-y)
\end{eqnarray*}

\begin{figure}[H]
\centering
\includegraphics*[scale=0.75]{sincos_nasaMetoda.png}
\caption{Površina funkcije $z=sin(x)+cos(x)$ interpolirane po naši metodi}
\label{fig:nasametoda}
\end{figure}

\section{Teoreti"cni del}
\begin{enumerate}
\item Najprej bomo dokazali, da ima funkcija $u(x,y)$ v robni to"cki (1,0) isto vrednost kot
funkcija $f_{10}$ v isti to"cki. Vzamemo funkcijo
\begin{eqnarray*}
    u(x,y) = f_{00}(x,y)q(x)q(y) + f_{10}(x,y)q(1-x)q(y) \\+ f_{01}(x,y)q(x)q(1-y) + f_{11}(x,y)q(1-x)q(1-y)
\end{eqnarray*}
in v njo vstavimo vrednosti $x=1$ ter $y=0$:
\begin{eqnarray*}
    u(1,0) = f_{00}(1,0)q(1)q(0) + f_{10}(1,0)q(0)q(0) \\+ f_{01}(x,y)q(1)q(1) + f_{11}(x,y)q(0)q(1)
\end{eqnarray*}

Polinom $q$ smo dolo"cili "ze v \ref{sec:utezna}.~poglavju, in sicer:
$$ q(x) = -6x^5 + 15x^4 - 10x^3 + 1 $$

$q(1) = 0$ ter $q(0) = 1$. Vidimo, da $f_{10}$ edini "clen, ki ni mno"zen s $q(1)$, torej nam na koncu ostane 

$$ u(1,0) = f_{10}(1,0) $$

Iz tega je razvidno, da imata $u(x,y)$ ter $f_{10}(x,y)$ v to"cki (1,0) iste vrednosti ter vrednosti odvodov.

\item 
V \ref{sec:utezna}.~poglavju smo pokazali, da je na"sa ute"zna funkcija polinom

$$ q(x) = -6x^5 + 15x^4 - 10x^3 + 1 $$.

Zdaj bomo dokazali, da je ute"z simetri"cna v smislu:

$$ q(1-x) = 1-q(x)$$.

Pokazali bomo, da je leva stran enaka desni strani. Desno stran bomo toliko
"casa spreminjali, da dobimo levo enako desni. Leva stran:

\begin{align*}
q(x) = &-6(1-x)^5 + 15(1-x)^4 -10(1-x)^3 + 1 = \\
     = &-6(-x^5 + 5x^4 - 10x^2 - 5x + 1) + \\ &15(x^4 - 4x^3 + 6x^2 - 4x + 1) - 10(x^2 -2x + 1) + 1 =\\
     = &6x^5 - 30x^3 + 60x^3 - 60x^2 + 30x - 6 + 15x^4 - 60x^3 + \\&90x^2 - 60x  + 15 + 10x^3 - 30x^2 + 30x - 10 + 1 =\\
     = &6x^5 - 15x^4 + 10x^3
\end{align*}

Desna stran:

\begin{align*}
1-q(x) &= 1-(-6x^5 + 15x^4 - 10x^3 + 1) = \\
       &= 1 + 6x^5 - 15x^4 + 10x^3 - 1 =\\
       &= 6x^5 - 15x^4 + 10x^3
\end{align*}
Vidimo, da sta leva in desna stran res enaki, zato trditev za to funkcijo iz danih pogojev velja.

To lahko v Mathematici poka"zemo v nekaj vrsticah:

\lstset{language=Mathematica}
\vspace{10pt}
\scriptsize
\begin{lstlisting}
k[x_] = a*x^5 + b*x^4 + c*x^3 + d*x^2 + e*x + f (*Definiramo polinom*)
Solve[{k[0] == 1, k[1] == 0, k'[0] == 0, k'[1] == 0, k''[0] == 0, 
    k''[1] == 0}, {a, b, c, d, e, f}](*Resimo sistem enacb*)
a = -6; b = 15; c = -10; d = 0; e = 0; f = 1; (*Spremenljivkam priredimo vrednosti*)
k[1 - x] == 1 - k[x] // Expand (*Testiramo enakost, test nam vrne True*)
\end{lstlisting}
\normalsize
\vspace{10pt}

Da sta funkciji res enaki lahko pogledamo še na graf. Rišemo z ukazom:
\vspace{10pt}
\scriptsize
\begin{lstlisting}
Plot[{k[1 - x], 1 - k[x], k[x]}, {x, -1, 2}, 
 PlotStyle -> {{Red, Thick}, {Yellow, Dashed, Thin}, {Blue, Thin}}]
\end{lstlisting}
\normalsize
\vspace{10pt}

\begin{figure}[H]
\centering
\includegraphics*[scale=0.65]{t.png}
\label{fig:sincos_lineinterp}
\end{figure}

Graf, na katerem imamo 3 funkcije. Rdeča funkcija je funkcija $q(1-x)$, rumena črtkana je $1-q(x)$, modra pa je osnovna funkcija $q(x)$. Slednjo preslikamo čez os $x$ in jo prestavimo za 1 višje, da dobimo prejšnji dve funkciji.

\item Dokazati moramo tudi, da naša utežna funkcija predstavlja razčlenitveno vsoto, ki je zaželjena lastnost, zato da bomo v primeru konstantnih podatkov dobili vedno dobili enako interpolacijsko funkcijo. V bistvu moramo dokazati trditev

$$ q(x)q(y) + q(1-x)q(y) + q(x)q(1-y) + q(1-x)q(1-y) = 1 $$

Ker smo "ze pokazali, da velja $q(1-x) = 1-q(x)$, lahko zgornjo ena"cbo raz"sirimo:

\begin{align*}
&q(x)q(y) + q(1-x)q(y) + q(x)q(1-y) + q(1-x)q(1-y) = \\
= &q(x)q(y) + (1-q(x))q(y) + q(x)(1-q(y)) + (1-q(x))(1-q(y)) = \\
= &q(x)q(y) + q(y) - q(x)q(y) + q(x) - q(x)q(y) + 1 - q(y) + q(x) + q(x)q(y) =\\
= &1
\end{align*}

Lahko testiramo "se v Mathematici:
\vspace{10pt}
\scriptsize
\begin{lstlisting}
k[x]*k[y] + k[1 - x]*k[y] + k[x]*k[1 - y] + k[1 - x]*k[1 - y] == 1 //Expand
(*Rezultat je True*)
\end{lstlisting}
\normalsize
\vspace{10pt}
\end{enumerate}

\appendix

\section{Primeri}

\begin{figure}[H]
\centering
\includegraphics*[scale=0.7]{randLinear.png}
\caption{Naključna funkcija A (linearna interpolacija)}
\label{fig:nasametoda}
\end{figure}

\begin{figure}[H]
\centering
\includegraphics*[scale=0.7]{randQOnly.png}
\caption{Naključna funkcija A (interpolacija samo z q(x))}
\label{fig:nasametoda}
\end{figure}

\begin{figure}[H]
\centering
\includegraphics*[scale=0.7]{randNasa.png}
\caption{Naključna funkcija A (naša metoda)}
\label{fig:nasametoda}
\end{figure}

\begin{figure}[H]
\centering
\includegraphics*[scale=0.65]{rand2Lin.png}
\caption{Naključna funkcija B (linearna interpolacija)}
\label{fig:nasametoda}
\end{figure}

\begin{figure}[H]
\centering
\includegraphics*[scale=0.65]{rand2Nasa.png}
\caption{Naključna funkcija B (naša metoda)}
\label{fig:nasametoda}
\end{figure}


\section{Koda}
\lstset{language=Matlab}
\vspace{10pt}
\scriptsize
\begin{lstlisting}
a=0; # da ne dojame kot funkcijo

# ----------------------------------------------------
# izracunaj utez
# ----------------------------------------------------
function [w]=utez(x)
  w=-6*x^5 + 15*x^4 - 10*x^3 + 0*x^2 + 0*x^1 + 1;
  #w=1-x;
endfunction

# ----------------------------------------------------
# ustvari matriko visin za GLAVNE TOCKE, ter matriki za X in Y 
# koordinate za vsako GLAVNO tocko,
# velikost je (stGlavnihKvY+3)*(stGlavnihKvX+3)
# ----------------------------------------------------
function [matrikaVisin,matrikaPozicijX,matrikaPozicijY]=napolniMatrikoGlavnihTock(
                                                         stGlavnihKvadratovX
                                                        ,minX
                                                        ,maxX)
  velikostGlavnegaKvX=(maxX-minX)/stGlavnihKvadratovX;
  stTockX=stGlavnihKvadratovX+1;
  matrika=zeros(stTockX,stTockX);
  offsetX=-(maxX-minX)/2;
  currentX=offsetX;
  currentY=offsetX;
  for i=[1:stTockX]
    for j=[1:stTockX]
      # tukaj lahko vnesemo poljubno funkcijo ce zelimo (zakomentiraj rand(..) spodaj)
      #matrikaVisin(i,j)=enacbaPovrsine(currentX,currentY);
      matrikaPozicijX(i,j)=currentX;
      matrikaPozicijY(i,j)=currentY;
      currentX+=velikostGlavnegaKvX;
    endfor
    currentX=offsetX;
    currentY+=velikostGlavnegaKvX;
  endfor  
  matrikaVisin=rand(stTockX,stTockX);
endfunction

# ----------------------------------------------------
# narisi mrezo
# ----------------------------------------------------
function narisi(mreza,stKvX,minX,maxX)
  tx = linspace (minX, maxX,stKvX+1);
  ty = linspace (minX, maxX,stKvX+1);
  [xx, yy] = meshgrid (tx, ty);
  mesh (tx, ty, mreza);
endfunction

# ----------------------------------------------------
# izracunaj odvod za dano tocko (za tocke glavnih kvadrator)
# ----------------------------------------------------
function [u_x,u_y,u_xx,u_yy,u_xy]=izracunajOdvodeZaDanoTocko(
                                     mrezaGlavnihKvadrator
                                    ,indexX
                                    ,indexY
                                    ,velGlKvX)
  sz=size(mrezaGlavnihKvadrator);  
  m=mrezaGlavnihKvadrator;
  h=velGlKvX;
  yIndexFac=1;
  
  if((indexX+1)>sz(1,2))
    indexX=sz(1,2)-1;
  endif
  
  if((indexX-1)<1)
    indexX=2;
  endif
  
  if((indexY+1)>sz(1,1))
    indexY=sz(1,1)-1;
  endif
  
  if((indexY-1)<1)
    indexY=2;
  endif
  
  u_x=m(indexY,indexX+1)-m(indexY,indexX-1);
  u_x=u_x/(2*h);
  
  u_y=m(indexY+1*yIndexFac,indexX)-m(indexY-1*yIndexFac,indexX);
  u_y=u_y/(2*h);
  
  u_xx=m(indexY,indexX+1)+m(indexY,indexX-1)-2*m(indexY,indexX);
  u_xx=u_xx/(h*h);
  
  u_yy=m(indexY+1*yIndexFac,indexX)+m(indexY-1*yIndexFac,indexX)-2*m(indexY,indexX);
  u_yy=u_yy/(h*h);
  
  u_xy=m(indexY+1*yIndexFac,indexX+1)-m(indexY+1*yIndexFac,indexX-1)- ...
       m(indexY-1*yIndexFac,indexX+1)+m(indexY-1*yIndexFac,indexX-1);
  u_xy=u_xy/(h*h);
  
endfunction

# ----------------------------------------------------
# funkcijske vrednosti za f00,f10,f01,f11
# normX in normY sta [0,1] x in y pa dejanske koordinate na robu
# ----------------------------------------------------
function [v]=f00(x,y,z,xv,yv,matrikaVisin,xi,yi,velKv)

  [dx,dy,ddx,ddy,dxy]=izracunajOdvodeZaDanoTocko(matrikaVisin,xi,yi,velKv);

  a=ddx/2;
  b=ddy/2;
  c=dxy;
  d=dx;
  e=dy;
  v=a*x^2+b*y^2+c*x*y+d*x+e*y+z;

endfunction
function [v]=f10(x,y,z,xv,yv,matrikaVisin,xi,yi,velKv)
  x=1-x;
  
  [dx,dy,ddx,ddy,dxy]=izracunajOdvodeZaDanoTocko(matrikaVisin,xi,yi,velKv);  
  
  a=ddx/2;
  b=ddy/2;
  c=-dxy;
  d=-dx;
  e=dy;
  v=a*x^2+b*y^2+c*x*y+d*x+e*y+z;
  
endfunction
function [v]=f01(x,y,z,xv,yv,matrikaVisin,xi,yi,velKv)
  y=1-y;
  
  [dx,dy,ddx,ddy,dxy]=izracunajOdvodeZaDanoTocko(matrikaVisin,xi,yi,velKv);

  a=ddx/2;
  b=ddy/2;
  c=-dxy;
  d=dx;
  e=-dy;
  v=a*x^2+b*y^2+c*x*y+d*x+e*y+z;
  
endfunction
function [v]=f11(x,y,z,xv,yv,matrikaVisin,xi,yi,velKv)
  x=1-x;
  y=1-y;
  
  [dx,dy,ddx,ddy,dxy]=izracunajOdvodeZaDanoTocko(matrikaVisin,xi,yi,velKv);

  a=ddx/2;
  b=ddy/2;
  c=dxy;
  d=-dx;
  e=-dy;
  v=a*x^2+b*y^2+c*x*y+d*x+e*y+z;
  
endfunction
# ----------------------------------------------------
# interpoliraj
# ----------------------------------------------------
function [matrikaVsehTock]=interpoliraj(matrikaVisin
                                        ,matrikaXCoord
                                        ,matrikaYCoord
                                        ,stGlKv
                                        ,stPodKv
                                        ,velikostGlKv)

  stVsehKv=stGlKv*stPodKv;
  matrikaVsehTock=zeros(stVsehKv+1,stVsehKv+1);
  size(matrikaVsehTock)
  # za vsak glavni kvadrat Y,X
  for gi=[1:stGlKv]
    for gj=[1:stGlKv]
    
      # iXY
      # | 00   01 |
      # | 10   11 |
      
      # indexi za v matriko glavnih tock (NE UPORABI ZA V PODTOCKE ZARADI ZAMIKA!)
      # giXX[x,y]
      gi00=[gj,gi];
      gi01=[gj,gi+1];
      gi10=[gj+1,gi];
      gi11=[gj+1,gi+1];
      
      x00=matrikaXCoord(gi00(1,2),gi00(1,1));  
      y00=matrikaYCoord(gi00(1,2),gi00(1,1));  
      x10=matrikaXCoord(gi10(1,2),gi10(1,1));  
      y10=matrikaYCoord(gi10(1,2),gi10(1,1));  
      x01=matrikaXCoord(gi01(1,2),gi01(1,1));  
      y01=matrikaYCoord(gi01(1,2),gi01(1,1));  
      x11=matrikaXCoord(gi11(1,2),gi11(1,1));  
      y11=matrikaYCoord(gi11(1,2),gi11(1,1));  
      
      v00=matrikaVisin(gi00(1,2),gi00(1,1));
      v01=matrikaVisin(gi01(1,2),gi01(1,1));
      v10=matrikaVisin(gi10(1,2),gi10(1,1));
      v11=matrikaVisin(gi11(1,2),gi11(1,1));
      
      
      
      stepSize=1/stPodKv;
      
      # normalized current values [0,1] 
      yParam=0;
      xParam=0;
      

      # za vsak podkvadrat Y,X
      for pi=[1:stPodKv+1]
        for pj=[1:stPodKv+1]
        

        
          # indexi za v matriko podtock (podkvadratov)
          # piXX[x,y]
          jj=gj-1;
          ii=gi-1;
          pi00=[jj*stPodKv+pj+0,ii*stPodKv+pi+0];
          pi01=[jj*stPodKv+pj+0,ii*stPodKv+pi+1];
          pi10=[jj*stPodKv+pj+1,ii*stPodKv+pi+0];
          pi11=[jj*stPodKv+pj+1,ii*stPodKv+pi+1];
                  
          w00x=utez(xParam);
          w00y=utez(yParam);
          
          w10x=utez(1-xParam);
          w10y=utez(yParam);

          w01x=utez(xParam);
          w01y=utez(1-yParam);
          
          w11x=utez(1-xParam);
          w11y=utez(1-yParam);
        
          
          interpX=(1-xParam)*x00+xParam*x10;
          interpY=(1-yParam)*y00+yParam*y01;
          

          u=f00(xParam,yParam,v00,x00,y00,matrikaVisin,
                gi00(1,1),gi00(1,2),velikostGlKv)*w00x*w00y+ ...
            f01(xParam,yParam,v01,x01,y01,matrikaVisin,
                gi01(1,1),gi01(1,2),velikostGlKv)*w01x*w01y+ ...
            f10(xParam,yParam,v10,x10,y10,matrikaVisin,gi10(1,1),
                gi10(1,2),velikostGlKv)*w10x*w10y+ ...
            f11(xParam,yParam,v11,x11,y11,matrikaVisin,gi11(1,1),
                gi11(1,2),velikostGlKv)*w11x*w11y;

          matrikaVsehTock(pi00(1,2),pi00(1,1))=u;

          xParam+=stepSize;
        endfor
        xParam=0;
        yParam+=stepSize;
      endfor
      
      
      
    endfor
  endfor
  
endfunction
# ----------------------------------------------------
# main
# ----------------------------------------------------

# MIN/MAX X IN Y MORATA BITI ENAKA, ST. X IN Y (POD)KV TUDI ENAKA

# ----------------------
#uporabnikovi parametri
# ----------------------
minimumX=-3; # zacetek X & Y
steviloKvadratovVSmeriX=4;
steviloPodkvadratorEnegaKvadrataVSmeriX=10; # 1 je isto kot "brez podkvadratov"


maximumX=-minimumX; # konec X
minimumY=minimumX; # zacetek Y
maximumY=maximumX; # konec Y
stVsehKvadratovVSmeriX=steviloKvadratovVSmeriX*steviloPodkvadratorEnegaKvadrataVSmeriX;

# ----------------------
# napolni matriko s vrednostmi visin, x coordinat in y coordinat
# (3 matrike, v vsaki ena komponenta (x,y,z)  za vsako tocko
# ----------------------

[matrikaVisinGlavnihTock,matrikaXCoord,matrikaYCoord]=napolniMatrikoGlavnihTock(
                                                    steviloKvadratovVSmeriX
                                                    ,minimumX
                                                    ,maximumX);

# ----------------------
# izracunaj interpolacijo
# ----------------------

velikostGlavnegaKv=(maximumX-minimumX)/steviloKvadratovVSmeriX;

matrikaVsehTock=interpoliraj(
                matrikaVisinGlavnihTock
                ,matrikaXCoord
                ,matrikaYCoord
                ,steviloKvadratovVSmeriX
                ,steviloPodkvadratorEnegaKvadrataVSmeriX
                ,velikostGlavnegaKv);


# ----------------------
# narisi
# ----------------------
#narisi(matrikaVisinGlavnihTock,steviloKvadratovVSmeriX,minimumX,maximumX);
narisi(matrikaVsehTock,stVsehKvadratovVSmeriX,minimumX,maximumX);
\end{lstlisting}
\normalsize
\vspace{10pt}
\end{document}
